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Abstract. A broad-beam-delivery system for heavy-charged-particle radiother- 
apy often employs multiple collimators and a range-compensating filter, which 
potentially offer complex beam customization. In treatment planning, it is how- 
ever difficult for a conventional pencil-beam algorithm to deal with these struc- 
tures due to beam-size growth during transport. This study aims to resolve the 
problem with a novel computational model. The pencil beams are initially de- 
fined at the range compensating filter with angular-acceptance correction for the 
upstream collimators followed by the range compensation effects. They are indi- 
vidually transported with possible splitting near the downstream collimator edges 
to deal with its fine structure. The dose distribution for a carbon-ion beam was 
calculated and compared with existing experimental data. The penumbra sizes 
of various collimator edges agreed between them to a submillimeter level. This 
beam-customization model will complete an accurate and efficient dose-calculation 
algorithm for treatment planning with heavy charged particles. 



PACS numbers: 87.53.Mr, 87.53.Pb, 87.53. Uv 
1. Introduction 

For heavy charged particle radiotherapy with protons and ions, broad-beam delivery 
methods (Coutrakon et al 1991, Kanai et al 1999) are mature technologies with 
persistent advantages of simplicity and robustness over emerging technologies of pencil- 
beam scanning methods (Lambert et al 2005). For a broad-beam system, a variety of 
volumetrically enlarged standard beams are prepared, among which an optimum one 
is applied to a given target. Target-specific customization is usually made with x-jaw, 
y-jaw, and multileaf collimators (XJC, YJC, and MLC) and custom-made accessories 
such as a patient collimator (PTC) and a range-compensating filter (RCF). While the 
downstream collimators form sharp field edges, the upstream collimators, which are 
mainly for radiation-protection purposes, form gentle field edges. Their combination 
will be useful for field patching techniques to form an irregular field with gently joining 
beams for improved robustness (Li et al 2007). 
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In treatment planning, a variety of pencil-beam (PB) algorithms are used for 
dose calculation (Hong et al 1996, Kanematsu et al 2006) despite intrinsic difficulty 
with the pencil beams that develop to overreach lateral heterogeneity (Goitein 1978, 
Petti 1992, Kohno et al 2004). For electron radiotherapy, the phase-space theory was 
rigorously applied to resolve the problem by periodical redefinition of ensemble of 
minimized pencil beams in the PB-redeffnition algorithm (Shiu and Hogstrom 1991). 
The principle of PB redefinition was applied to heavy charged particles to address the 
effects of multiple collimators in the monochromatic PB approximation (Kanematsu 
et al 2008b). However, its rigorous application to a heterogeneous system requires 
polychromatic energy spectra, which would be computationally demanding for heavy 
charged particles with sharp Bragg peaks. It will be thus difficult to cope with range 
compensation or patient heterogeneity in that approach. 

Recently, Kanematsu et al (2009) proposed an alternative approach, the PB- 
splitting algorithm, where monochromatic pencil beams dynamically split into smaller 
ones near a lateral density interface. Automatically, fine pencil beams are densely 
arranged only where they are necessary while otherwise large pencil beams are sparsely 
arranged for efficient dose calculation. In conjunction with the grid-dose-sprcading 
convolution (Kanematsu et al 2008a), the PB-splitting algorithm demonstrated 
feasibility of accurate patient dose calculation while minimizing the impact of recursive 
beam multiplication (Kanematsu 2011). 

In this study, we further extend the PB-splitting approach to beam-customization 
devices to deal with their physical structures accurately and efficiently and to complete 
a consistent algorithmic framework for dose calculation in treatment planning. In 
the following sections, we define the model elements that were mostly diverted from 
previous studies, construct a novel and original beam-customization model, and 
examine its validity for a test-beam experiment. 

2. Materials and methods 

2.1. Model elements 

2.1.1. Beam source and beam's eye view A beam source is defined as the best 
approximate point from which radiating particles will have the same fluence reduction 
with distance. The formulation differs among beam-spreading methods and often 
between the transverse x and y axes, i.e., at height z$ x for x and zs y for y. The 
particles incoming to a point in the field, which is normally the isocenter, are 
projected back onto the x and y source planes to define rms source sizes as x and 
(j gy in the Gaussian approximation. Although a range-modulated beam should be 
ideally subdivided into energy components of different source heights and sizes, it is 
approximately represented by a single component of average behavior in this study. 

A beam's-eye-view (BEV) image is defined as an n x m matrix of S x S square- 
sized pixels starting at (x\,yi) on the isocenter plane. For BEV pixel ij, pixel position 
(xj,yi) and the line connecting to the x and y sources are defined as 




(1) 
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2.1.2. Collimator Following the thick-collimator model (Kanematsu et al 2006), two 
identical apertures on the top and bottom faces are associated with every collimator, 
which are modeled as two-dimensional bitmaps. Matrix T a describes aperture a with 
n a x m a elements of transmission T aia j a = 1 (transmit) or (block) for i a e [l,n a ] 
and j a G [1 , m a ] . The pixel-i a j a position is given by 

(Xaj a ,VaiJ = {Xal,Val) + (ja ~ Ma ~ l)<*a, (2) 

where (x al ,y al ) and <5 a are the first pixel position and the square pixel size of the 
bitmap image. For arbitrary point (x a ,y a ) on the aperture plane, intersecting pixel 
i a j a is determined with the nearest integer function |_ ] as 

Da Va\ 



Ja — 



+ 1, ia = 



Sa 



+ 1, (3) 



Sa 

and the distance to the nearest aperture edge, 

d aia3a = t .,.,, rr {SaV(i'a-ia) 2 + (j' a - ja) 2 ) , (4) 

is quickly referenced from the distance map filled by the distance-transform algorithm 
(Borgefors 1986). 

2.1.3. Range compensating filter A RCF made of a tissue-like material of effective 
density pc is similarly described by an nc x mc matrix of range shifts Sq, first pixel 
position (xci,yci)t & n d pixel size Sc- In this study, we deal with a single RCF of a 
flat downstream face at height zq. The stopping and scattering effects of the RCF are 
approximated by a local interaction at the midpoint of the beam path in the structure 
(Gottschalk et al 1993), i.e., at height z — zq + 0.5 Sa c j c /pc for RCF pixel icjc- 

2.1. 4- Pencil beam Following the original PB-splitting algorithm (Kanematsu et 
al 2009), the present PB model is based on the Fermi-Eyges theory (Eyges 1948) 
for stopping and scattering (Kanematsu 2009, Gottschalk 2010) excluding hard 
interactions that are implicitly included in the depth-dose curve. 

A Gaussian pencil beam is characterized by position r, direction v, number of 
particles n, residual range R, and phase-space variances of the projected angle 9 and 
transverse displacement t, which develop in a tissue-like medium by step As as 

Ar = vAs, Av = 0, An = 0, Ai?=-pAs, (5) 

= 1%-0.16 ( M ~ - 92 ln R (6) 

1000 y \m v ) R + AR' y ' 



AO I -- \ e 2 '+ ^Ad 2 } As. Al' 2 = 



■2 01 + i d 2 '+^A9 2 } As 



As, (7) 



where p is the stopping-power ratio of the medium to water (Kanematsu et al 2003) 
and m/m p and q are the particle mass and charge in units of those of a proton. 

To limit excessive beam multiplication, pencil beams subject to splitting should 
have sufficient particles, i.e., n/n > K n , where tiq is the number for the original beam 
and cutoff. When a pencil beam of rms size cr t = V(t 2 ) spreads beyond the 

lateral density interface at distance di nt from the beam center, it splits into M x M 
daughter beams downsized by factor om as 

2 for cr 2 K d a t < d int < n d a t ( n 2 \ ( V3/2 

M = { 3 for a 3 n d a t < d [nt <a 2 K d a tl \ a 3 \ = \ l/y/2 \ , (8) 

4 for dint < o- 3 K d a t 
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where K<j is a parameter that limits the fraction of overreaching particles. With respect 
to the mother beam, daughter a/3 (a, j3 € [1, M]) is downscaled, displaced, redirected, 
and downsized while conserving focal distance t 2 /8t and local mean square angle 
2 - M 2 /¥, as 

' 1/2 1/2 
1/4 1/2 1/4 
1/8 3/8 3/8 1/8 

-l 



n a /3 — fhlafMfj ' 



Ar, 



a/3, 




t 2 a/3 — 0~M t 2 , 




9t aP = a 2 M 6t, 9 2 a0 = 9 2 - (1 - a 2 M ) 0t /t 2 , 



(9) 
(10) 

(11) 
(12) 



where fM a fMp is the share of daughter a/3 and [iM a and \Im$ are the displacement 
factors for transverse directions e t rs and e u ~ e v . 



2.2. The beam- delivery system model 

2.2.1. Pencil-beam generation For every BEV pixel ij, pencil beam b is placed and 
the PB parameters are defined at the effective interaction point in the RCF as 

2 2 1 _1 / 2 

2-) +1 

ZSx ) \ZSyJ 



~0i 



z c + 0.5Sa c j c /p< 



i'b 



V, 



1 



rib = n 0b = $Qij $ 



R = Rq, 



(13) 
(14) 
(15) 
(16) 

(17) 

y/ZSx — ZOijy/ ZSy — Zoij Z Sx Z Sy 11 

where RCF intersection pixel zcjc is determined in analogy with ([1]) and (j3]) and 
open-field fluence $or/ and range Rq are usually given by measurement. 



w 1 

2 \ z Sx - z 



et 
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ZSy - Zq 1 



ZQij ZSy — ZQij (S 2 



2.2.2. Upstream collimation The upstream collimators restrict angular acceptance 
of particles incoming to each PB origin. A pencil beam will be fully blocked when it 
is far away from any one of the apertures or fully transmitted when it is in the middle 
of all of them, or 



T h = 



for { 3 a | (z a > Z Q ); (T aiaja = 0) A [d aiaJa > 3cr gb (z a - z )] } 

1 for { y a \ (z a > z );(T aiaja = l) A [d aiaJa > 3o- gb (z a - z )]}, 



(18) 



where T b is the transmission factor of beam 6, i a j a is the beam intersection pixel 
of aperture a, and factor 3 to rms projected angle crg b = \/(9 2 b) secures three 
standard deviations for edge distance d a i a j a ■ For partial transmission, we calculate the 
geometrical acceptance of particles incoming to the PB origin. As shown in figure [TJ 




Figure 1. Definition of geometrical parameters (angle 9 and axes x, z, s, and t) 
in the z— x view, where symbols S, b, and p indicate the source, a pencil beam, 
and a particle, respectively. 



with small orthogonal angles and <fi about the PB axis, constituent-particle direction 
v p is defined as 

Vp = (vpx? v p y i v p^) ~ v b + tan 6< e x + tan0 e y , (19) 
which translates into geometrical line (6,<fi), or 

| x{z) = x + {z- z ) v Px /v Pz ^ 
\ y{z) = 2/o + (z - z ) v Py /v Pz ' 

Only particles passing through all the apertures can get to the PB origin to redefine 
the number of particles, the direction, and the mean square angle as 

n b = n 0b T bl T b = 11 ^ 2 e " \.l T *i'J' a > ( 21 ) 

b {a\z a >z } 



v bz ' v bz J T b JJ_ n 2na g 2 



b {a\z a >z a } 

Bh -T b L,^?: 11^-2-. ( 5 

where is the aperture pixel in which line (0, <j>) intersects. In practice, these 
integrals are made numerically at 0.2 ag b sampling intervals for ±3crg b regions. 



2.2.3. Range compensation The RCF shortens the residual range of the pencil beam 
by the thickness of the intersecting pixel as AR b = —Sa c j c and increases the mean 
square angle by A8 2 b in (|6]) before the beam is transported downstream. 



2.2.4- Downstream collimation Every pencil beam is individually transported by 
(|5 |) — (|7 [) through downstream apertures. At an aperture, which is practically either 
the top or bottom face of an optional PTC, pencil beams near the edge will be 
partially transmitted. Incidentally, edge distance d aia j a in (0| naturally corresponds to 
density-interface distance d nlt in ([8]) for PB splitting. At every downstream aperture, 
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Figure 2. Beam-customization devices of the experiment; (a) the y = cross- 
section view, (b) the beam's eye view on the isocenter plane, where the filled 
areas represent XJC, YJC, MLC, and PTC from upstream to downstream and 
the hatched area represents RCF. 



multiplicity M is appropriately determined while limiting overreaching particles to 
below 2% by setting n n — 3. In the case of splitting, daughter beams are defined 
according to (j9|)- (|12|) and then individually transported downstream starting from 
the current aperture with possible recursive splitting in the same manner. The pencil 
beams that are finally out of the aperture will be blocked by setting rib = 0, which 
addresses the partial-blocking effect of the collimator. 

2.3. Experimental validation 

2.3.1. Apparatus As this study shares the objective of beam-customization modeling 
with the former study in the PB-redefinition approach (Kanematsu et al 2008b), we 
use the same experimental data, where a broad carbon-ion beam of residual range 
Ro = 19.6 cm in water was customized with an XJC at height 117-137 cm, a YJC 
at 96-116 cm, and a partially effective MLC at 69-83 cm, a 3-cm PMMA half-plate 
RCF at 35-38 cm, and an 8-cm-square PTC at 22-27 cm as shown in figure [21 Four 
lines of the in-air dose profiles on the isocenter plane were measured along the x axis 
at y = —2 cm and 1 cm and along the y axis at x = — 1 cm and 1 cm. The 20%-80% 
penumbra sizes (« 1.68 at) were 0.58 cm for the XJC edge and 0.48 cm for the YJC 
edge, which translate into rms source sizes asx = (0.58/1.68)(940 — 1 17)/1 17 = 2.43 
cm at z Sx = 940 cm and a Sy = (0.48/1. 68)(1040 - 96)/96 = 2.81 cm at z Sy = 1040 
cm. The tissue-air ratio for the 3-cm PMMA (p = 1.16) was measured to be 0.951. 

2.3.2. Implementation In the calculation, 100 x 100 dose grids in a single layer were 
arranged on the isocenter plane at 1-mm intervals. The open field of uniform fluence 
($o = 1) was subdivided into the BEV image pixels of size S = 0.5 mm on the 
isocenter plane, to each of which a pencil beam was defined at the effective scattering 
point of the RCF. For every pencil beam, upstream collimation by the XJC, the YJC, 
and the MLC, range shift and scattering by the RCF, and beam transport including 
collimation and splitting by the PTC down to the isocenter plane were applied. The 
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Figure 3. Calculated dose distribution in gray scale for the customized carbon- 
ion field, where the dotted contours represent 20%, 50%, and 80% dose levels and 
the gray lines indicate the dose-profiling positions. 



in-air dose distribution on the isocenter plane was calculated with 

Tl, D (x-x b ) 2 + (i;-H b ) 2 ( 1 f n 

"" % fi ^1 . D* b = H 

^ 0.951 



D(x,y) 



h l b 



19.6 cm 



for Ri, 16.1 cm, 



(24) 



where D^ b is the tissue-air ratio or the dose per fluence for beam b. 

To verify the effectiveness of PB splitting for the PTC edge, we calculated dose 
distributions at heights cm (isocenter plane) and 20 cm (immediate downstream) 
by relocating the dose grids and compared them with corresponding non-splitting 
calculations, for which we disabled splitting by setting Kd = 0. 



3. Results 



In the calculation, 40000 beams were originally defined at the RCF, 23912 of them 
passed through the upstream collimators, and 20444 of them passed through the PTC 
to end up with 36704 dose-contributing beams by splitting. The CPU time of a 2.4 
GHz Intel Core 2 Duo processor amounted to 1.30 s and 1.25 s for the calculations 
with and without PB splitting. Figure [3] shows the calculated dose distribution. The 
dip and bump along the y axis are attributed to scattering by the PMMA half plate. 
Sharpness of the field edge was strongly correlated with the distance to the effective 
collimator. 

In the experiment, the uncertainty of the scanned detector positions was 0.1 mm 
and that of the collimator positions was < 0.5 mm according to the specifications. 
The latter may only shift the edge position and will not influence the penumbra size. 
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Figure 4. Dose profiles along the x axis at y = —2 cm (a), 3/ = +1 cm (b), along 
the y axis at x = —2 cm (c), and x = +1 cm (d), where the solid lines are the 
calculations and the open circles are the measurements. 



The single-point dose uncertainty was evaluated to be 0.3% in repeated measurements, 
which is negligible for penumbra analysis. Figure HI shows the calculated and measured 
doses profiles, where the measured doses are in fact the dose ratios of the customized 
field to the open field to compensate for the fluence non-uniformity. Unexpectedly, the 
customized-field doses were higher than the open-field doses by a few percent. That 
may be attributed to the contribution of particles hard-scattered by the collimators, 
which was not considered in the present model. 

From these profiles, the 20%-80% penumbra sizes were obtained by reading 20% 
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Table 1. Measured and calculated 20%-80% penumbra sizes for the customized 
carbon-ion field. 
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Figure 5. Dose profiles for PTC edge; (a) immediate downstream with RCF, (b) 
immediate downstream without RCF, (c) isocenter plane with RCF, (d) isocenter 
plane without RCF, where the solid and the dashed lines are calculations with 
and without PB splitting, respectively, and the open circles are measurements. 



and 80% dose positions by linear interpolation of two sampling points, which brings 
dominant uncertainty amounting to a fraction of the sampling interval of 1 mm. The 
measured penumbra sizes were then corrected to quadratically exclude 1.68cr s i ze with 
effective dosimeter size <7 s i ze = 0.5 mm for a 2-mm<ji> pinpoint chamber. Table [1] 
summarizes the resultant penumbra sizes. These measurements and calculations 
agreed to a submillimeter level, which is consistent with the estimated uncertainty. 

The effectiveness of beam splitting for the PTC-edge sharpening is shown in 
figure [SJ where panels (c) and (d) show enlarged views of panels (c) and (d) in figure [4] 
with additional lines for non-splitting calculations. The PB splitting reasonably 
sharpened the field edges at the immediate downstream and made better agreement 
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with the measurements on the isocenter plane. Ironically, the contamination of 
collimator-scattered particles happened to compensate substantially for the lack of 
edge sharpening in the tail regions. 

4. Discussion 

It is one of the algorithmic novelties of this study to originate the pencil beams at the 
effective scattering points of the RCF regardless of upstream collimation. Then, the 
lateral heterogeneity of the RCF is naturally irrelevant to the minimized pencil beams. 
Upstream collimation is reasonably modeled as filtering of particles in the angular 
distribution to correct the phase-space parameters of the defined pencil beams. 

While the PB size and density are generally arbitrary in PB algorithms, small 
size and high density are required to represent sharp edges of downstream collimation. 
In the present model, the sharp PTC edge was naturally realized by splitting of the 
pencil beams. In the former study (Kanematsu et al 2008b), because the pencil 
beams could not be redefined as monochromatic after range compensation, they 
were only artificially downsized for edge sharpening. The downsizing strength was 
empirically determined to reproduce the 20%-80% penumbra size on the isocenter 
plane while overlooking the other aspects. In fact, while the resultant penumbra 
sizes were equivalcntly good for both models, the dose profiles in the upstream was 
unphysically bouncy in the former study due to insufficient density, which could be 
clinically problematic. 

In the original PB-splitting algorithm (Kanematsu et al 2009), the overreaching 
condition was defined as the one-standard-deviation distance (k^ = 1) to a 10% density 
change. That was because its objective heterogeneity was moderate density variation 
among body tissues. This study deals with solid and precisely defined collimator 
edges, for which the «d = 3 distance to an aperture edge may be more appropriate. 

In the present example, the PTC was effective only for approximately 1/4 of the 
field edge. The PB splitting was limited to the pencil beams around the effective 
edge and actually increased the number of beams by 80% and the CPU time by 
4%. This discrepancy is mainly attributed to computational overhead for generation, 
upstream collimation, and range compensation of the pencil beams. Although we 
only dealt with the planer grids in this case, the PB splitting would not add severe 
computational load even for volumetric grids when used with the grid-dose-spreading 
convolution (Kanematsu 2011). 

In heavy charged particle radiotherapy, target doses are predominantly formed 
by Bragg peaks of primary particles. Hard-scattered particles are generally out of 
the scope of practical PB algorithms due to difficulty in their modeling. Fortunately, 
the collimator-scattered particles tend to lose large energy in the collimator and thus 
naturally attenuate with depth (van Luijk et al 2001). Nevertheless, Kimstrand et al 
(2008) included the collimator-scatter contribution in a convolution algorithm using 
Monte-Carlo-generated kernels. That approach may be valid and will further improve 
the accuracy if combined with the present model. 

5. Conclusions 

We have developed a calculation model for customization of a broad beam of heavy 
charged particles based on the PB-splitting algorithm. In this model, a broad beam is 
decomposed into pencil beams of various size that is necessarily and sufficiently small 
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to deal with structures of the beam-customization devices accurately and efficiently. 
Also, placement of the PB origins at the effective scattering points in the RCF 
effectively reduced the relevant heterogeneity and greatly simplified the algorithm 
using only monochromatic pencil beams. 

The performance of the model was tested against existing experimental data, 
which demonstrated that the penumbra size for various collimator edges in a single field 
was accurate to a submillimeter level. This beam-customization part can be naturally 
combined with the patient-dose-calculation part that is similarly based on the PB- 
splitting algorithm (Kanematsu 2011) to complete an accurate and efficient dose 
calculation algorithm for treatment planning of heavy-charged-particle radiotherapy. 
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